



cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
      subroutine k11(o_n,A,na,B,nb,c)
      implicit none
      integer na,nb
Cf2py intent(in) na,nb
      real*8 c
Cf2py intent(in) c
      real*8 A(na),B(nb)
Cf2py intent(in) A,B
Cf2py depend(na) A
Cf2py depend(nb) B
      real*8 o_n(na*nb)
Cf2py intent(out) o_n
Cf2py depend(na,nb) o_n
      integer i,j
      o_n=0
      do i=1,na,1
      do j=1,nb,1
      o_n(j+(i-1)*nb)=c*A(i)*B(j)
      end do
      end do
c      return
      end subroutine k11
c
c
      subroutine k21r(o_n,A,na,B,nbx,nby,c)
      implicit none
      integer na,nbx,nby
Cf2py intent(in) na,nbx,nby
      real*8 A(na),B(1:nbx,1:nby)
Cf2py intent(in) A,B
Cf2py depend(na) A
Cf2py depend(nbx,nby) B
      real*8 c
Cf2py intent(in) c
      real*8 o_n(1:nbx,1:na*nby)
Cf2py intent(out) o_n
Cf2py depend(na,nbx,nby) o_n
      integer i,i2,j2
      o_n=0
      do i=1,na,1
      do i2=1,nbx,1
      do j2=1,nby,1
      o_n(i2,j2+(i-1)*nby)=A(i)*B(i2,j2)*c
      end do
      end do
      end do
c      return
      end subroutine k21r
c
c
      subroutine k21c(o_n,A,B,c,na,nbx,nby)
      implicit none
      integer, intent(in) :: na,nbx,nby
Cf2py intent(in) na,nbx,nby
      real*8, intent(out), dimension(1:na*nbx,1:nby) :: o_n
Cf2py intent(out) o_n
Cf2py depend(na,nbx,nby) o_n
      real*8, intent(in) :: c
Cf2py intent(in) c

      real*8 A(na),B(1:nbx,1:nby)
Cf2py intent(in) A,B
Cf2py depend(na) A
Cf2py depend(nbx,nby) B
      integer i,i2,j2
      o_n=0
      do i=1,na,1
      do i2=1,nbx,1
      do j2=1,nby,1
      o_n(i2+(i-1)*nbx,j2)=A(i)*B(i2,j2)*c
      end do
      end do
      end do
c      return
      end subroutine k21c
c
c
      subroutine k22(o_n,A,nax,nay,da,B,nbx,nby,db,c)
      implicit none
      integer nax,nay,nbx,nby,da,db
Cf2py intent(in) nax,nay,nbx,nby
      real*8 A(1:nax,1:nay),B(1:nbx,1:nby)
Cf2py intent(in) A,B
Cf2py depend(nax,nay) A
Cf2py depend(nbx,nby) B
      real*8 c
Cf2py intent(in) c
      real*8 o_n(1:nax*nbx,1:nay*nby)
Cf2py intent(out) o_n
Cf2py depend(nax,nay,nbx,nby) o_n
      integer i1,i2,j1,j2,io,jo
      o_n=0

      do i1=1,nax,1
        if (da.eq.0) then
          do j1=1,nay,1
      do i2=1,nbx,1
        if (db.eq.0) then
          do j2=1,nby,1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end do
        else if (db.eq.1) then
          j2=i2
          io=i2+(i1-1)*nbx
          jo=j2+(j1-1)*nby
          o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
        else if (db.eq.2) then
          if (i2.eq.1) then
            j2=2
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else if (i2.eq.nbx) then
            j2=nbx-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else
            j2=i2+1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
            j2=i2-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end if
        end if
      end do
          end do
        else if (da.eq.1) then
          j1=i1
      do i2=1,nbx,1
        if (db.eq.0) then
          do j2=1,nby,1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end do
        else if (db.eq.1) then
          j2=i2
          io=i2+(i1-1)*nbx
          jo=j2+(j1-1)*nby
          o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
        else if (db.eq.2) then
          if (i2.eq.1) then
            j2=2
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else if (i2.eq.nbx) then
            j2=nbx-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else
            j2=i2+1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
            j2=i2-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end if
        end if
      end do
        else if (da.eq.2) then
          if (i1.eq.1) then
            j1=2
      do i2=1,nbx,1
        if (db.eq.0) then
          do j2=1,nby,1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end do
        else if (db.eq.1) then
          j2=i2
          io=i2+(i1-1)*nbx
          jo=j2+(j1-1)*nby
          o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
        else if (db.eq.2) then
          if (i2.eq.1) then
            j2=2
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else if (i2.eq.nbx) then
            j2=nbx-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else
            j2=i2+1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
            j2=i2-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end if
        end if
      end do
          else if (i1.eq.nax) then
            j1=nax-1
      do i2=1,nbx,1
        if (db.eq.0) then
          do j2=1,nby,1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end do
        else if (db.eq.1) then
          j2=i2
          io=i2+(i1-1)*nbx
          jo=j2+(j1-1)*nby
          o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
        else if (db.eq.2) then
          if (i2.eq.1) then
            j2=2
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else if (i2.eq.nbx) then
            j2=nbx-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else
            j2=i2+1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
            j2=i2-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end if
        end if
      end do
          else
            j1=i1+1
      do i2=1,nbx,1
        if (db.eq.0) then
          do j2=1,nby,1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end do
        else if (db.eq.1) then
          j2=i2
          io=i2+(i1-1)*nbx
          jo=j2+(j1-1)*nby
          o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
        else if (db.eq.2) then
          if (i2.eq.1) then
            j2=2
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else if (i2.eq.nbx) then
            j2=nbx-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else
            j2=i2+1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
            j2=i2-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end if
        end if
      end do
            j1=i1-1
      do i2=1,nbx,1
        if (db.eq.0) then
          do j2=1,nby,1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end do
        else if (db.eq.1) then
          j2=i2
          io=i2+(i1-1)*nbx
          jo=j2+(j1-1)*nby
          o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
        else if (db.eq.2) then
          if (i2.eq.1) then
            j2=2
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else if (i2.eq.nbx) then
            j2=nbx-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else
            j2=i2+1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
            j2=i2-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end if
        end if
      end do
          end if
        end if
      end do

c      return
      end subroutine k22
c
c
c
c
c

      subroutine kk22(o_n,na,nb,nc,da,db,dc,A,B,C,g)
c kron(a,kron(b,c)
      implicit none
      integer na,nb,nc,da,db,dc
Cf2py intent(in) na,nb,nc,da,db,dc
      real*8 A(1:na,1:na),B(1:nb,1:nb),C(1:nc,1:nc)
Cf2py intent(in) A,B,C
Cf2py depend(na) A
Cf2py depend(nb) B
Cf2py depend(nc) C
      real*8 g
Cf2py intent(in) g
      real*8 o_n(1:na*nb*nc,1:na*nb*nc)
Cf2py intent(out) o_n
Cf2py depend(na,nb,nc) o_n
      integer ia,ib,ic,ja,jb,jc,io,jo
      o_n=0
c      do ia=1,na,1
c      do ib=1,nb,1
c      do ic=1,nc,1
c      do ja=1,na,1
c      do jb=1,nb,1
c      do jc=1,nc,1
c      io=ic+(ib-1)*nc+(ia-1)*nb*nc
c      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
c      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
c      end do
c      end do
c      end do
c      end do
c      end do
c      end do

      do ia=1,na,1
        if (da.eq.0) then
          do ja=1,na,1
      do ib=1,nb,1
        if (db.eq.0) then
          do jb=1,nb,1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end do
        else if (db.eq.1) then
          jb=ib

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
        else if (db.eq.2) then
          if (ib.eq.1) then
            jb=2

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else if (ib.eq.nb) then
            jb=nb-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else
            jb=ib+1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
            jb=ib-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end if
        end if
      end do
          end do
        else if (da.eq.1) then
          ja=ia
      do ib=1,nb,1
        if (db.eq.0) then
          do jb=1,nb,1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end do
        else if (db.eq.1) then
          jb=ib

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
        else if (db.eq.2) then
          if (ib.eq.1) then
            jb=2

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else if (ib.eq.nb) then
            jb=nb-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else
            jb=ib+1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
            jb=ib-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end if
        end if
      end do
        else if (da.eq.2) then
          if (ia.eq.1) then
            ja=2
      do ib=1,nb,1
        if (db.eq.0) then
          do jb=1,nb,1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end do
        else if (db.eq.1) then
          jb=ib

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
        else if (db.eq.2) then
          if (ib.eq.1) then
            jb=2

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else if (ib.eq.nb) then
            jb=nb-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else
            jb=ib+1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
            jb=ib-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end if
        end if
      end do
          else if (ia.eq.na) then
            ja=na-1
      do ib=1,nb,1
        if (db.eq.0) then
          do jb=1,nb,1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end do
        else if (db.eq.1) then
          jb=ib

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
        else if (db.eq.2) then
          if (ib.eq.1) then
            jb=2

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else if (ib.eq.nb) then
            jb=nb-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else
            jb=ib+1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
            jb=ib-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end if
        end if
      end do
          else
            ja=ia+1
      do ib=1,nb,1
        if (db.eq.0) then
          do jb=1,nb,1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end do
        else if (db.eq.1) then
          jb=ib

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
        else if (db.eq.2) then
          if (ib.eq.1) then
            jb=2

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else if (ib.eq.nb) then
            jb=nb-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else
            jb=ib+1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
            jb=ib-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end if
        end if
      end do
            ja=ia-1
      do ib=1,nb,1
        if (db.eq.0) then
          do jb=1,nb,1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end do
        else if (db.eq.1) then
          jb=ib

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
        else if (db.eq.2) then
          if (ib.eq.1) then
            jb=2

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else if (ib.eq.nb) then
            jb=nb-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else
            jb=ib+1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
            jb=ib-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end if
        end if
      end do
          end if
        end if
      end do
c      return
      end subroutine kk22











c
c
      subroutine sk22(o_n,na,nb,ca,cb,da,db,dia,dib,A,B,IA,IB,g)
c o=kron22(A,IB)
c o+=kron22(IA,B)
      implicit none
      integer na,nb,da,db,dia,dib
Cf2py intent(in) na,nb
      real*8 A(1:na,1:na),B(1:nb,1:nb)
      real*8 IA(1:na,1:na),IB(1:nb,1:nb)
Cf2py intent(in) A,B,IA,IB
Cf2py depend(na) A,IA
Cf2py depend(nb) B,IB
      real*8 ca,cb,g
Cf2py intent(in) ca,cb,g
      real*8 o_inter(1:na*nb,1:na*nb), o_n(1:na*nb,1:na*nb)
Cf2py intent(out) o_n
Cf2py depend(na,nb) o_n,o_inter
      o_n=0
      o_inter=0
      call k22(o_inter,A,na,na,da,IB,nb,nb,dib,ca)
      call k22(o_n,IA,na,na,dia,B,nb,nb,db,cb)
      o_n=(o_n+o_inter)*g
c      return
      end subroutine sk22
c
c
      subroutine sskk22(o_n,na,nb,nc,ca,cb,cc,d,A,B,C,IA,IB,IC,g)
c o=kk22(A,kk22(IB,IC))
c o+=kk22(IA,kk22(B,IC))
c +=kk22(IA,kk22(IB,C))
c kron(a,kron(b,c)
      implicit none
      integer na,nb,nc,d(6)
Cf2py intent(in) na,nb,nc,d
      real*8 A(1:na,1:na),B(1:nb,1:nb),C(1:nc,1:nc)
      real*8 IA(1:na,1:na),IB(1:nb,1:nb),IC(1:nc,1:nc)
Cf2py intent(in) A,B,C,IA,IB,IC
Cf2py depend(na) A,IA
Cf2py depend(nb) B,IB
Cf2py depend(nc) C,IC
      real*8 ca,cb,cc,g
Cf2py intent(in) ca,cb,cc,g
      real*8 o_n(1:na*nb*nc,1:na*nb*nc)
      real*8 o_inter1(1:na*nb*nc,1:na*nb*nc)
      real*8 o_inter2(1:na*nb*nc,1:na*nb*nc)
Cf2py intent(out) o_n
Cf2py depend(na,nb,nc) o_n
      o_n=0
      o_inter1=0
      o_inter2=0
      call kk22(o_inter1,na,nb,nc,d(1),d(5),d(6),A,IB,IC,ca)
      call kk22(o_inter2,na,nb,nc,d(4),d(2),d(6),IA,B,IC,cb)
      call kk22(o_n,na,nb,nc,d(4),d(5),d(3),IA,IB,C,cc)
      o_n=g*(o_n+o_inter1+o_inter2)
c      return
      end subroutine sskk22
c
c
      subroutine sp(o_n,na,nb,cx,cy,cz,d,Ax,Ay,Az,Bx,By,Bz,c)
      implicit none
c scalar product vector of n component of tensor rank da db
c declaration of dimensions
      integer na,nb,d(6)
Cf2py intent(in) na,nb,d
      real*8 c,cx,cy,cz
Cf2py intent(in) c,cx,cy,cz
      real*8 o_n(1:na*nb,1:na*nb)
      real*8 o_inter1(1:na*nb,1:na*nb),o_inter2(1:na*nb,1:na*nb)
Cf2py intent(out) o_n
Cf2py depend(na,nb) o_n
      real*8 Ax(1:na,1:na),Ay(1:na,1:na),Az(1:na,1:na)
Cf2py intent(in) Ax,Ay,Az
Cf2py depend(na) Ax,Ay,Az
      real*8 Bx(1:nb,1:nb),By(1:nb,1:nb),Bz(1:nb,1:nb)
Cf2py intent(in) Bx,By,Bz
Cf2py depend(nb) Bx,By,Bz
      o_inter1=0
      o_inter2=0
      o_n=0
      call k22(o_inter1,Ax,na,na,d(1),Bx,nb,nb,d(4),cx)
      call k22(o_inter2,Ay,na,na,d(2),By,nb,nb,d(5),cy)
      call k22(o_n,Az,na,na,d(3),Bz,nb,nb,d(6),cz)
      o_n=c*(o_n+o_inter1+o_inter2)
c      return
      end subroutine sp

      subroutine sptv(o_n,na,Ax,Ay,Az,B,c)
      implicit none
c scalar product vector of n component of tensor rank da db
c declaration of dimensions
      integer na
Cf2py intent(in) da
      real*8 o_n(1:na,1:na)
Cf2py intent(out) o_n
Cf2py depend(na) o_n
      real*8 Ax(1:na,1:na),Ay(1:na,1:na),Az(1:na,1:na)
Cf2py intent(in) Ax,Ay,Az
Cf2py depend(na) Ax,Ay,Az
      real*8 B(3)
      real*8 c
Cf2py intent(in) B,c
      o_n=0
      o_n=B(1)*Ax
      o_n=o_n+B(2)*Ay
      o_n=o_n+B(3)*Az
      o_n=o_n*c
c      return
      end subroutine sptv



cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
      subroutine build_sx(o,n)
      implicit none
c a subroutine to build a Sz matrix with a spin n=2s+1
      integer n,k
Cf2py intent(in) n
      real*8 o(1:n,1:n)
      real*8 Ip(1:n,1:n)
Cf2py intent(out) o
Cf2py depend(n) o,Ip
      real*8 ms
      real*8 lm(n)

      Ip = 0
      o=0

      lm=0
      ms=-(n-1)/2.0
      do k=1,n,1
      lm(k) = ms
      ms=ms+1
      end do

      do k=1,n-1,1
      Ip(k , k + 1) = sqrt(((n-1)/2.0 - lm(k)) * ((n-1)/2.0 + lm(k)+ 1))
      end do

      o = (1.0 / 2) * (Ip + transpose(Ip))
c      return
      end subroutine build_sx
c
c
c
      subroutine build_sy(o,n)
      implicit none
c a subroutine to build a Sz matrix with a spin n=2s+1
      integer n,k
Cf2py intent(in) n
      real*8 o(1:n,1:n)
      real*8 Ip(1:n,1:n)
Cf2py intent(out) o
Cf2py depend(n) o,Ip
      real*8 ms
      real*8 lm(n)

      Ip = 0
      o=0

      lm=0
      ms=-(n-1)/2.0
      do k=1,n,1
      lm(k) = ms
      ms=ms+1
      end do

      do k=1,n-1,1
      Ip(k , k + 1) = sqrt(((n-1)/2.0 - lm(k)) * ((n-1)/2.0 + lm(k)+ 1))
      end do

      o = (1.0 / 2) * (-Ip + transpose(Ip))
c      return
      end subroutine build_sy
c
c
      subroutine build_sz(o,n)
      implicit none
c a subroutine to build a Sz matrix with a spin n=2s+1
      integer n,k
Cf2py intent(in) n
      real*8 o(1:n,1:n)
Cf2py intent(out) o
Cf2py depend(n) o
      real*8 ms

      o=0

      ms=(n-1)/2.0
      do k=1,n,1
      o(k,k)=ms
      ms=ms-1
      end do
c      return
      end subroutine build_sz
c
c
      subroutine build_sd(o,n)
      implicit none
c a subroutine to build an identity matrix with a dimension n
      integer n,k
Cf2py intent(in) n
      real*8 o(1:n,1:n)
Cf2py intent(out) o
Cf2py depend(n) o
      o=0
      do k=1,n,1
      o(k,k)=1
      end do
c      return
      end subroutine build_sd
c
c
      subroutine bfield(o,b0,xyz,f,phi,A,t)
      real*8 o(3)
Cf2py intent(out) o
      real*8 xyz(3)
Cf2py intent(in) xyz
      real*8 b0,f,phi,A,t
Cf2py intent(in) b0,f,phi,A,t
      real*8 pi,norm

      pi=3.1415026535
      o=0
      norm=sqrt(xyz(1)**2+xyz(2)**2+xyz(3)**2)
      o = (b0 + A* sin(2*pi*t*f+phi)) / norm * xyz
c      return
      end subroutine bfield
c
c
      function lorentzian(i,c,w,x)
      implicit none
c w the width
c c the center
      real*8 lorentzian
      real*8 i,c,w,x
Cf2py intent(in) i,c,w,x
      real*8 num,denom,pi
      num=0
      denom=0
      pi=3.1415926535
      num=i * (w / (2 * pi))
      denom=((x-c)*(x-c)+(w/2)*(w/2))
      lorentzian=num/denom
      end function lorentzian
c
c
      function gaussian(i,c,w,x)
      implicit none
c w the width
c c the center
      real*8 gaussian,pi
      real*8 i,c,w,x
Cf2py intent(in) i,c,w,x
      pi=3.1415926535
      gaussian=i/(w*sqrt(2*pi))*exp(-(c-x)*(c-x)/(2*w*w))
      end function gaussian

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine build_ci(o,d2,d3,d4,d5,Cp,Cap)
      implicit none
c DECLARATIONS of VARIABLES
c declaration of dimensions
      integer d2,d3,d4,d5
Cf2py intent(in) d2,d3,d4,d5
c
c declaration of the matrix
c  output matrix
      real*8 o(1:d2,1:d2)
Cf2py intent(out) o
Cf2py depend(d2) o
c  capture part of the o matrix
      real*8 Cm11(1:d3*d4,1:d3*d4),Cm21(1:d3*d4,1:d3*d4)
      real*8 Cm12(1:d2,1:d2),Cm22(1:d2,1:d2)
      real*8 Sz(1:2,1:2),Sd(1:2,1:2),Id(1:d5,1:d5)
c
c declaration of some constant
      real*8 Cp,Cap,Calpha,Cbeta
Cf2py intent(in) Cp,Cap
c
c DEFINITIONS and INITIALIZATIONS
c  base rates
      Calpha = (Cp - Cap) / 2.0
      Cbeta = (Cp + Cap) / 2.0
c matrix
      Sz=0
      Sd=0
      Id=0
      Cm11=0
      Cm12=0
      Cm21=0
      Cm22=0
      call build_sz(Sz,2)
      call build_sd(Sd,2)
      call build_sd(Id,d5)
      call k22(Cm11,Sz,2,2,1,Sz,2,2,1,1)
      call k22(Cm12,Cm11,4,4,1,Id,d5,d5,1,1)
      call k22(Cm21,Sd,2,2,1,Sd,2,2,1)
      call k22(Cm22,Cm21,4,4,1,Id,d5,d5,1)

      o = Calpha *  Cm12 + Cbeta * Cm22

c      return
      end subroutine build_ci
c
c
      subroutine build_cj(o,Mbj2,d2,d3,d4,d5,Cp,Cap)
      implicit none
c DECLARATIONS of VARIABLES
c declaration of dimensions
      integer d2,d3,d4,d5
Cf2py intent(in) d2,d3,d4,d5
c
c declaration of the matrix
c  output matrix
      real*8 o(1:d2,1:d2)
Cf2py intent(out) o
Cf2py depend(d2) o
c  coefficient matrix (element wise squared)
      real*8 Mbj2(1:d2,1:d2)
Cf2py intent(in) Mbj2
Cf2py depend(d2) Mbj2
c  capture part of the o matrix
      real*8 Cv1(d3*d4),Cv2(d5),Cv3(d2)
c
c declaration of some constant
      real*8 Cp,Cap
Cf2py intent(in) Cp,Cap
c
c declaration of counters
      integer i,j
c DEFINITIONS and INITIALIZATIONS
c matrix
      Cv2=1
      Cv1(1)=Cp
      Cv1(2)=Cap
      Cv1(3)=Cap
      Cv1(4)=Cp
      Cv3=0
      call k11(Cv3,Cv1,d3*d4,Cv2,d5,1)
      do i=1,d2,1
      do j=1,d2,1
      o(i,j)=Cv3(j)*Mbj2(j,i)
      end do
      end do
c      return
      end subroutine build_cj
c
c
      subroutine build_rj(o,Mbj2,d2,d3,d4,d5,Rp,Rap)
      implicit none
c DECLARATIONS of VARIABLES
c declaration of dimensions
      integer d2,d3,d4,d5
Cf2py intent(in) d2,d3,d4,d5
c
c declaration of the matrix
c  output matrix
      real*8 o(1:d2,1:d2)
Cf2py intent(out) o
Cf2py depend(d2) o
c  coefficient matrix (element wise squared)
      real*8 Mbj2(1:d2,1:d2)
Cf2py intent(in) Mbj2
Cf2py depend(d2) Mbj2
c  capture part of the o matrix
      real*8 Rv1(d3*d4),Rv2(d5)
      real*8 Rv3(d2)
Cf2py depend(d3,d4) Rv1
Cf2py depend(d5) Rv2
Cf2py depend(d2) Rv3
c
c declaration of some constant
      real*8 Rp,Rap
Cf2py intent(in) Rp,Rap
c
c declaration of counters
      integer i,j
c DEFINITIONS and INITIALIZATIONS
c matrix
      Rv1(1)=Rp
      Rv1(2)=Rap
      Rv1(3)=Rap
      Rv1(4)=Rp
      Rv2=1
      Rv3=0
      call k11(Rv3,Rv1,d3*d4,Rv2,d5,1)
      do i=1,d2,1
      do j=1,d2,1
      o(i,i)=o(i,i)+Rv3(j)*Mbj2(j,i)
c      print*,o(i,i)
      end do
      end do
c      print*,o
c      return
      end subroutine build_rj
c
c
      subroutine build_rk(o,Mbj2,d2,d3,d4,d5,Rp,Rap)
      implicit none
c DECLARATIONS of VARIABLES
c declaration of dimensions
      integer d2,d3,d4,d5
Cf2py intent(in) d2,d3,d4,d5
c
c declaration of the matrix
c  output matrix
      real*8 o(1:d5,1:d2)
Cf2py intent(out) o
Cf2py depend(d2,d5) o
c  coefficient matrix (element wise squared)
      real*8 Mbj2(1:d2,1:d2)
Cf2py intent(in) Mbj2
Cf2py depend(d2) Mbj2
c  capture part of the o matrix
      real*8 Rv1(d3*d4),Rv2(d5)
      real*8 Rv3(d2)
c
c declaration of some constant
      real*8 Rp,Rap
Cf2py intent(in) Rp,Rap
c
c declaration of counters
      integer i,j,k
c DEFINITIONS and INITIALIZATIONS
c matrix
      o=0
      Rv1(1)=Rp
      Rv1(2)=Rap
      Rv1(3)=Rap
      Rv1(4)=Rp
      Rv3=0

      do i=1,d5,1
      Rv2=0
      Rv2(i)=1
      call k11(Rv3,Rv1,d3*d4,Rv2,d5,1)
      do j=1,d2,1
      do k=1,d2,1
      o(i,j)=o(i,j)+Rv3(k)*Mbj2(k,j)
      end do
      end do
      end do
c      return
      end subroutine build_rk
c
c
      subroutine build_pek(o,d5,PE)
      implicit none
c DECLARATIONS of VARIABLES
c declaration of dimensions
      integer d5
Cf2py intent(in) d5
c
c declaration of the matrix
c  output matrix
      real*8 o(1:d5,1:d5)
Cf2py intent(out) o
Cf2py depend(d5) o
c  coefficient matrix (element wise squared)
c
c declaration of some constant
      real*8 PE
Cf2py intent(in) PE
c
c declaration of counters
      integer i
c DEFINITIONS and INITIALIZATIONS
c matrix
      o=0
      do i=1,d5,1
      o(i,i)=PE / d5
      end do
c      return
      end subroutine build_pek
c
c
      subroutine build_pei(o,d2,d3,d4,d5,PE)
      implicit none
c DECLARATIONS of VARIABLES
c declaration of dimensions
      integer d2,d3,d4,d5
Cf2py intent(in) d1,d2,d3,d4,d5
c
c declaration of the matrix
c  output matrix
      real*8 o(1:d2,1:d5)
Cf2py intent(out) o
Cf2py depend(d2,d5) o
c  photo excitation part of the o matrix
      real*8 PEv1(d3*d4)
      real*8 PEm1(1:d2,1:d5)
      real*8 Id(1:d5,1:d5)
c declaration of some constant
      real*8 PE
Cf2py intent(in) PE
c
      integer i,j
c
c
c matrix
      o=0
      PEv1=1
      PEm1=0
      Id=0
      call build_sd(Id,d5)
c      call k21c(PEm1,PEv1,d3*d4,Id,d5,d5,1)
      do i=1,d2,1
      do j=1,d5,1
      o(i,j)=PEm1(i,j)*PE/d2
      end do
      end do
c      print *,o
c      return
      end subroutine build_pei
c
c
      subroutine nonresmat(o,Mbj2,d1,d2,d3,d4,d5,Cp,Cap,Rp,Rap,PE)
      implicit none
c DECLARATIONS of VARIABLES
c declaration of dimensions
      integer d1,d2,d3,d4,d5
Cf2py intent(in) d1,d2,d3,d4,d5
c declaration of boundaries for the matrix
      integer c1i,c1e,c2i,c2e,c3i,c3e
c
c declaration of the matrix
c  output matrix
      real*8 o(1:d1,1:d1)
Cf2py intent(out) o
Cf2py depend(d1) o
c  coefficient matrix (element wise squared)
      real*8 Mbj2(1:d2,1:d2)
Cf2py intent(in) Mbj2
Cf2py depend(d2) Mbj2

c Matrix related to Ci\Cj
      real*8 Ci(1:d2,1:d2),Cj(1:d2,1:d2)
Cf2py depend(d2) Ci,Cj
      real*8 Cp,Cap,Calpha,Cbeta
Cf2py intent(in) Cp,Cap
      real*8 Cm11(1:d3*d4,1:d3*d4),Cm21(1:d3*d4,1:d3*d4)
      real*8 Cm12(1:d2,1:d2),Cm22(1:d2,1:d2)
      real*8 Cv1(d3*d4),Cv2(d5),Cv3(d2)
Cf2py depend(d3,d4) Cm11,Cm21,Cv1
Cf2py depend(d2) Cm12,Cm22,Cv3
Cf2py depend(d5) Cv2

c Matrix related to Rj\Rk
      real*8 Rj(1:d2,1:d2),Rk(1:d5, 1:d2)
Cf2py depend(d2) Rj
Cf2py depend(d2,d5) Rk
      real*8 Rv1(d3*d4),Rv2(d5)
      real*8 Rv3(d2)
Cf2py depend(d3,d4) Rv1
Cf2py depend(d5) Rv2
Cf2py depend(d2) Rv3
      real*8 Rp,Rap
Cf2py intent(in) Rp,Rap

c Matrix related to PEi\PEk
      real*8 PEi(1:d2,1:d5),PEk(1:d5,1:d5)
Cf2py depend(d2,d5) PEi
Cf2py depend(d5) PEk
      real*8 PE
Cf2py intent(in) PE
      real*8 PEv1(1:d3*d4)
      real*8 PEm1(1:d2,1:d5)

c Operator matrixw
      real*8 Sz(1:2,1:2),Sd(1:2,1:2),Id(1:d5,1:d5)

      integer i,i2,io,j,j2,k,l,m

c DEFINITIONS and INITIALIZATIONS
c constants
c  dimensions
      c1i=1
      c1e=d2
      c2i=c1e+1
      c2e=c1e+d2
      c3i=c2e+1
      c3e=c2e+d5

c matrix
      o=0
      Ci=0
      Cj=0
      Rj=0
      Rk=0
      PEi=0
      PEk=0
      Sz=0
      Sd=0
      Id=0

      call build_sz(Sz,2)
      call build_sd(Sd,2)
      call build_sd(Id,d5)

c Construct self.NRrates
c Ci: assuming there is no coupling at this step
c(even if we write the transformation matrix)
cccccccccccccccccccccccccccccccccccccccccccc
c
c      call build_ci(Ci,d2,d3,d4,d5,Cp,Cap)
c
cccccccccccccccccccccccccccccccccccccccccccc

      Calpha = (Cp - Cap) / 2.0
      Cbeta = (Cp + Cap) / 2.0
c matrix
      Cm11=0
      Cm12=0
      Cm21=0
      Cm22=0

c
c
      o(c1i:c1e, c1i:c1e) = -Ci
c      o(c2i:c2e, c1i:c1e) = Cj
c      o(c2i:c2e, c2i:c2e) = -Rj
c      o(c3i:c3e, c2i:c2e) = Rk
c      o(c3i:c3e, c3i:c3e) = -PEk
c      o(c1i:c1e, c3i:c3e) = PEi
      o=real(o)

c      return
      end subroutine nonresmat
c
c
      subroutine resmat(o,dtm,f,w,mwe,ev,d1,d2)
c
      implicit none
c DECLARATIONS of VARIABLES
c declaration of dimensions
      integer d1,d2
Cf2py intent(in) d1,d2
c
c declaration of the matrix
c  output matrix
      real*8 o(1:d1,1:d1)
Cf2py intent(in,out) o
Cf2py depend(d1) o
c  coefficient matrix (element wise squared)
      real*8 dtm(1:d2,1:d2)
Cf2py intent(in) dtm
Cf2py depend(d2) dtm
      real*8 ev(d2)
Cf2py intent(in) ev
Cf2py depend(d2) ev
c
c declaration of some constant
      real*8 f,w,mwe
Cf2py intent(in) f,w,mwe
      real*8 dts,gaussian,intensity,center
c dipolar_transition_strength
      integer i,j,im,jm


c
      do i=1,d2,1
      do j=1,d2,1
      center = abs(ev(i) - ev(j))
      intensity = dtm(i,j)**2
c the dipolar transition matrix element has to be squared to get intensity
      dts= gaussian(intensity,center,w,f)
      im=d2+i
      jm=d2+j
c      if (center.lt.9082) then
c      if (center.gt.9070) then
c      print*,dts, intensity,center, w, f
c      end if
c      end if

      o(im,im)=o(im,im)-mwe*dts
      o(jm,jm)=o(jm,jm)-mwe*dts
      o(im,jm)=o(im,jm)+mwe*dts
      o(jm,im)=o(jm,im)+mwe*dts
      dts=0
      intensity=0
      end do
      end do
      o=real(o)
c      return
      end subroutine resmat


      subroutine test(PEi,d2,d3,d4,d5,PE)
      implicit none
c DECLARATIONS of VARIABLES
c declaration of dimensions
      integer d2,d3,d4,d5
Cf2py intent(in) d2,d3,d4,d5
c declaration of boundaries for the matrix
c declaration of the matrix
c  output matrix
c Matrix related to PEi\PEk
      real*8 PEi(1:d2,1:d5)
Cf2py intent(out) PEi
Cf2py depend(d2,d5) PEi
      real*8 PE
Cf2py intent(in) PE
      real*8 PEv1(1:d3*d4)
      real*8 PEm1(1:d2,1:d5)

c Operator matrixw
      real*8 Id(1:d5,1:d5)

      integer i,i2,j,j2

      external :: build_sd
      external :: k21c

c DEFINITIONS and INITIALIZATIONS
c constants
c  dimensions

c matrix
      PEi=0
      Id=0
      call build_sd(Id,d5)
cccccccccccccccccccccccccccccccccccccccccccc
c
c      call build_pei(PEi,d2,d3,d4,d5,PE)
c
cccccccccccccccccccccccccccccccccccccccccccc
      PEv1=1
      PEm1=0
      print*,PEm1
cccccccccccccccc
      call k21c(PEm1,PEv1,Id,1,d3*d4,d5,d5)
c      do i=1,d3*d4,1
c      do i2=1,d5,1
c      do j2=1,d5,1
c      PEm1(i2+(i-1)*d5,j2)=PEv1(i)*Id(i2,j2)
c      end do
c      end do
c      end do
c      end call k21c(PEm1,PEv1,Id,1,d3*d4,d5,d5)
cccccccccccccccc
      print*,PEm1
      do i=1,d2,1
      do j=1,d5,1
      PEi(i,j)=PEm1(i,j)*PE/d2
      end do
      end do

c      return
      end subroutine test
